

This example implements the Cahn-\/\-Hilliard equation for phase-\/field modeling of a two species, as described by the following weak form of the P\-D\-E. The scalar fields are composition one, $c_1$, and composition two, $c_2$. Note the application of the higher-\/order Dirichlet boundary conditions $\nabla c_i\cdot\boldsymbol{n}=0$ using Nitsche's method.

Cahn-\/\-Hilliard\-:

\begin{eqnarray*} 0 &=& \int_\Omega \left(w_1\frac{c_1 - c_{1_{prev}}}{\mathrm{d}t} + M\left(\nabla w_1\cdot(f_{,c_1c_1}\nabla c_1) + f_{,c_1c_2}\nabla c_2) + \kappa_1\nabla^2 w_1\nabla^2 c_1\right)\right) dV\\ &\phantom{=}& - \int_{\partial\Omega} \left(w_1j_n + M\kappa_1\left(\nabla^2c_1(\nabla w_1\cdot\boldsymbol{n}) + \nabla^2w_1(\nabla c_1\cdot\boldsymbol{n})\right) - \tau(\nabla w_1\cdot\boldsymbol{n})(\nabla c_1\cdot\boldsymbol{n})\right) dS\\ 0 &=& \int_\Omega \left(w_2\frac{c_2 - c_{2_{prev}}}{\mathrm{d}t} + M\left(\nabla w_2\cdot(f_{,c_2c_1}\nabla c_1) + f_{,c_2c_2}\nabla c_2) + \kappa_1\nabla^2 w_2\nabla^2 c_2\right)\right) dV\\ &\phantom{=}& - \int_{\partial\Omega} \left(w_2j_n + M\kappa_1\left(\nabla^2c_2(\nabla w_2\cdot\boldsymbol{n}) + \nabla^2w_2(\nabla c_2\cdot\boldsymbol{n})\right) - \tau(\nabla w_2\cdot\boldsymbol{n})(\nabla c_2\cdot\boldsymbol{n})\right) dS \end{eqnarray*}

Free energy density\-:

\begin{eqnarray*} f(c_1,c_2) = c_2(c_1-0.1)^2(c_1-0.9)^2 + (c_2-0.2)^2(c_2-0.8)^2 \end{eqnarray*}

\section*{Implementation\-: Level 1 users }

To implement this model, we will specify the following through defining user functions\-: \par

\begin{DoxyItemize}
\item Initial conditions \par

\item Constitutive model (via free energy density functions) \par

\item Parameter values \par

\item Weak form of the P\-D\-Es \par

\end{DoxyItemize}

First, we include the header file declaring the required user functions. These functions will be defined in this file.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


Now, we first define any optional user functions. Optional user functions have a default definition that can be redefined by the user using a function pointer. This will be done in the {\ttfamily define\-Parameters} function. The available list of optional user functions includes\-: {\ttfamily boundary\-Conditions}, {\ttfamily scalar\-Initial\-Conditions}, {\ttfamily vector\-Initial\-Conditions}, {\ttfamily load\-Step}, {\ttfamily adaptive\-Time\-Step}, and {\ttfamily project\-Fields}. In this example, we redefine only the {\ttfamily scalar\-Initial\-Conditions} function, while using the default functions for the others.

{\bfseries  The {\ttfamily scalar\-Initial\-Conditions} function }

We initialized the two composition fields to be random about 0.\-5.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


{\bfseries  Free energy density derivative functions }

This phase-\/field implementation requires the first and second derivatives of the chemical free energy density function $f(c_1,c_2) = c_2(c_1-0.1)^2(c_1-0.9)^2 + (c_2-0.2)^2(c_2-0.8)^2$. We define the function computing $\partial^2 f/\partial c \partial c$ here. Note that this free energy derivative function is used only in this file. It is not a member of any class, nor will we use it to set any function pointers.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


{\bfseries  The {\ttfamily define\-Parameters} function }

The user is required to define the {\ttfamily define\-Parameters} and {\ttfamily residual} functions. The {\ttfamily define\-Parameters} defines variables and functions in the {\ttfamily App\-Ctx} object. The {\ttfamily App\-Ctx} object is defined in the app\-Ctx.\-h file. This function is used to define any values in {\ttfamily user} that will be needed in the problem. It is also used to set any function pointers for user functions that we have redefined.

Many of these values can be overwritten by the parameters.\-prm file, which we will look at later.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


Here, we define the mesh by setting the number of elements in each direction, e.\-g. a 100x100 element mesh.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


We also define the dimensions of the domain, e.\-g. a unit square.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


We define the initial time step and total simulation time. We also have the options to use restart files, in which case we would set the iteration index and time at which to start. We leave these values at zero to begin a new simulation. We also have the option to output results at regular intervals (e.\-g. every 5 time steps).


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


We specify the number of vector and scalar solution and projection fields by adding the name of each field to their respective vector. Here, we have one scalar solution field (the composition). We do not use any vector solution fields or projection fields in this example.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


We can specify the polynomial order of the basis splines, as well as the global continuity. Note that the global continuity must be less than the polynomial order. Here, we use quadratic basis functions with C-\/1 global continuity.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


Finally, we redirect the desired user function pointers to the {\ttfamily scalar\-Initial\-Conditions} function that we defined above. This completes the {\ttfamily define\-Parameters} function.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


{\bfseries  The {\ttfamily residual} function }

The residual function defines the residual that is to be driven to zero. This is the central function of the code. It is set up to follow the analytical weak form of the P\-D\-E. It has a number of arguments that give problem information at the current quadrature point.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


{\ttfamily d\-V} is a boolean, \char`\"{}true\char`\"{} if {\ttfamily residual} is being called for the volume integral and \char`\"{}false\char`\"{} if {\ttfamily residual} is being called for the surface integral.\par
{\ttfamily d\-S} is a boolean, \char`\"{}false\char`\"{} if {\ttfamily residual} is being called for the volume integral and \char`\"{}true\char`\"{} if {\ttfamily residual} is being called for the surface integral.\par
{\ttfamily x} gives the coordinates of the quadrature point.\par
{\ttfamily normal} gives the unit normal for a surface quadrature point.\par
{\ttfamily c} gives the information (values, gradients, etc.) for the scalar solution fields at the current quadrature point (see documentation for solution\-Scalars class).\par
{\ttfamily u} gives the information (values, gradients, etc.) for the vector solution fields at the current quadrature point (see documentation for solution\-Vectors class).\par
{\ttfamily w1} gives the information for the scalar test functions.\par
{\ttfamily w2} gives the information for the vector test functions.\par
{\ttfamily user} is a structure available for parameters related to the initial boundary value problem (e.\-g. elasticity tensor).\par
{\ttfamily r} stores the scalar value of the residual for the weak form of the P\-D\-E which is then used by the core assembly functions.

The following functions are available for the solution objects {\ttfamily c} and {\ttfamily u}, where the argument is the field index, i.

{\ttfamily c.\-val(i)} -\/ Value of scalar field i, scalar \par
{\ttfamily c.\-grad(i)} -\/ Gradient of scalar field i, 1st order tensor \par
{\ttfamily c.\-hess(i)} -\/ Hessian of scalar field i, 2nd order tensor \par
{\ttfamily c.\-laplacian(i)} -\/ Laplacian of scalar field i, scalar \par
{\ttfamily c.\-val\-P(i)} -\/ Value of scalar field i at previous time step, scalar \par
{\ttfamily c.\-grad\-P(i)} -\/ Gradient of scalar field i at previous time step, 1st order tensor \par
{\ttfamily c.\-hess\-P(i)} -\/ Hessian of scalar field i at previous time step, 2nd order tensor \par
{\ttfamily c.\-laplacian\-P(i)} -\/ Laplacian of scalar field i at previous time step, scalar

{\ttfamily u.\-val(i)} -\/ Value of vector field i, 1st order tensor \par
{\ttfamily u.\-grad(i)} -\/ Gradient of vector field i, 2nd order tensor \par
{\ttfamily u.\-hess(i)} -\/ Hessian of vector field i, 3rd order tensor \par
{\ttfamily u.\-val\-P(i)} -\/ Value of vector field i at previous time step, 1st order tensor \par
{\ttfamily u.\-grad\-P(i)} -\/ Gradient of vector field i at previous time step, 2nd order tensor \par
{\ttfamily u.\-hess\-P(i)} -\/ Hessian of vector field i at previous time step, 3rd order tensor

Similar functions are available for the test functions. Also, the following tensor operations are useful\-:

Tensor operations\-: \par
{\ttfamily operator+} -\/ tensor addition \par
{\ttfamily operator-\/} -\/ tensor subraction \par
{\ttfamily operator$\ast$} -\/ single contraction between tensors or scalar multiplication \par
{\ttfamily double\-\_\-contract} -\/ double contraction of two 2nd order tensors, or a 4th order tensor and a 2nd order tensor. \par
{\ttfamily trans( )} -\/ transpose 2nd order tensor \par
{\ttfamily trace( )} -\/ trace of 2nd order tensor \par
{\ttfamily det( )} -\/ determinant of 2nd order tensor \par
{\ttfamily inv( )} -\/ inverse of 2nd order tensor \par


The example code here implements the weak form for the Cahn-\/\-Hilliard equations, as shown above.

First, we set the values for necessary parameters, using some predefined material parameters.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


Next, we get the values for the free energy derivatives based on the current quadrature point.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


Now, we compute the residual in a manner very similar to the analytical form\-:

\begin{eqnarray*} 0 &=& \int_\Omega \left(w_1\frac{c_1 - c_{1_{prev}}}{\mathrm{d}t} + M\left(\nabla w_1\cdot(f_{,c_1c_1}\nabla c_1) + f_{,c_1c_2}\nabla c_2) + \kappa_1\nabla^2 w_1\nabla^2 c_1\right)\right) dV\\ &\phantom{=}& - \int_{\partial\Omega} \left(w_1j_n + M\kappa_1\left(\nabla^2c_1(\nabla w_1\cdot\boldsymbol{n}) + \nabla^2w_1(\nabla c_1\cdot\boldsymbol{n})\right) - \tau(\nabla w_1\cdot\boldsymbol{n})(\nabla c_1\cdot\boldsymbol{n})\right) dS\\ 0 &=& \int_\Omega \left(w_2\frac{c_2 - c_{2_{prev}}}{\mathrm{d}t} + M\left(\nabla w_2\cdot(f_{,c_2c_1}\nabla c_1) + f_{,c_2c_2}\nabla c_2) + \kappa_1\nabla^2 w_2\nabla^2 c_2\right)\right) dV\\ &\phantom{=}& - \int_{\partial\Omega} \left(w_2j_n + M\kappa_1\left(\nabla^2c_2(\nabla w_2\cdot\boldsymbol{n}) + \nabla^2w_2(\nabla c_2\cdot\boldsymbol{n})\right) - \tau(\nabla w_2\cdot\boldsymbol{n})(\nabla c_2\cdot\boldsymbol{n})\right) dS \end{eqnarray*}


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


Finally, we include a file that instatiates the template functions {\ttfamily define\-Parameters} and {\ttfamily residual}. This bit of code will generally be the same for any problem (unless you decide to use a different automatic differentation library); the user does not need to modify it.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


The complete implementation can be found at \href{https://github.com/mechanoChem/mechanoChemIGA/blob/master/initBounValProbs/CahnHilliard_twoSpecies/2D/userFunctions.cc}{\tt Github}.

\section*{Parameters file\-: Interface for level 2 users }

Now let's look at the parameters file, {\ttfamily parameters.\-prm}. The advantages of the parameters file are that these values can be changed without recompiling the code and it can provide a clean interface to the code. 

The parameters defined in the parameters file overwrite any previous values defined in the {\ttfamily define\-Parameters} function. Anything following the pound sign (\#) is a comment. A parameter is defined using the syntax\-:

{\ttfamily set} {\ttfamily parameter\-Name} {\ttfamily =} {\ttfamily parameter\-Value} 

There is a set list of variables that can be read from the parameters file. Anything else will be added to the {\ttfamily mat\-Param} structure with a double number type. Tensor objects can follow the format\-: 1 x 1 or \mbox{[}1,1\mbox{]} or (1,1), where the number of components must equal the spatial dimension of the problem.

In this example file, we begin by specifying the spatial dimension, the geometry dimensions, and the mesh size\-:


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


We then define time stepping, restart information, output frequency, and spline parameters.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


Note that we don't need to include all (or even any) of these parameters in this file. We defined default values previously.

The complete parameters file can be found at \href{https://github.com/mechanoChem/mechanoChemIGA/blob/master/initBounValProbs/CahnHilliard_twoSpecies/2D/parameters.prm}{\tt Github}. 